To be able to edit code and run cells, you need to run the notebook yourself. Where would you like to run the notebook?

In the cloud (experimental)

Binder is a free, open source service that runs scientific notebooks in the cloud! It will take a while, usually 2-7 minutes to get a session.

On your computer

(Recommended if you want to store your changes.)

  1. Copy the notebook URL:
  2. Run Pluto

    (Also see: How to install Julia and Pluto)

  3. Paste URL in the Open box

Frontmatter

If you are publishing this notebook on the web, you can set the parameters below to provide HTML metadata. This is useful for search engines and social media.

Author 1

homework 9, version 0

👀 Reading hidden code
226 μs

Submission by: Jazzy Doe (jazz@mit.edu)

👀 Reading hidden code
30.7 ms

Homework 9: Climate modeling I

18.S191, fall 2020

👀 Reading hidden code
296 μs
👀 Reading hidden code
# edit the code below to set your name and kerberos ID (i.e. email without @mit.edu)

student = (name = "Jazzy Doe", kerberos_id = "jazz")

# you might need to wait until all other cells in this notebook have completed running.
# scroll around the page to see what's up
15.5 μs

Let's create a package environment:

👀 Reading hidden code
230 μs
Package dependencies
👀 Reading hidden code
begin
import Pkg
Pkg.activate(mktempdir())
Pkg.add([
"Plots",
"PlutoUI",
"LaTeXStrings",
"Distributions",
"Random",
])
using LaTeXStrings
using Plots
using PlutoUI
using Random, Distributions
Random.seed!(123)
md"##### Package dependencies"
end
❔
  Activating new project at `/tmp/jl_4MXcLy`
    Updating registry at `~/.julia/registries/General.toml`
   Resolving package versions...
    Updating `/tmp/jl_4MXcLy/Project.toml`
  [31c24e10] + Distributions v0.25.118
  [b964fa9f] + LaTeXStrings v1.4.0
  [91a5bcdd] + Plots v1.40.14
  [7f904dfe] + PlutoUI v0.7.64
  [9a3f8284] + Random
    Updating `/tmp/jl_4MXcLy/Manifest.toml`
  [6e696c72] + AbstractPlutoDingetjes v1.3.2
  [66dad0bd] + AliasTables v1.1.3
  [d1d4a3ce] + BitFlags v0.1.9
  [d360d2e6] + ChainRulesCore v1.25.1
  [9e997f8a] + ChangesOfVariables v0.1.10
  [944b1d66] + CodecZlib v0.7.8
  [35d6a980] + ColorSchemes v3.29.0
  [3da002f7] + ColorTypes v0.12.1
  [c3611d14] + ColorVectorSpace v0.11.0
  [5ae59095] + Colors v0.13.1
  [34da2185] + Compat v4.16.0
  [f0e56b4a] + ConcurrentUtilities v2.5.0
  [187b0558] + ConstructionBase v1.5.8
  [d38c429a] + Contour v0.6.3
  [9a962f9c] + DataAPI v1.16.0
  [864edb3b] + DataStructures v0.18.22
  [b429d917] + DensityInterface v0.4.0
  [31c24e10] + Distributions v0.25.118
  [ffbed154] + DocStringExtensions v0.9.5
  [460bff9d] + ExceptionUnwrapping v0.1.11
  [c87230d0] + FFMPEG v0.4.2
  [1a297f60] + FillArrays v1.13.0
  [53c48c17] + FixedPointNumbers v0.8.5
  [1fa38f19] + Format v1.3.7
  [28b8d3ca] + GR v0.73.6
  [42e2da0e] + Grisu v1.0.2
  [cd3eb016] + HTTP v1.10.16
  [34004b35] + HypergeometricFunctions v0.3.28
  [47d2ed2b] + Hyperscript v0.0.5
  [ac1192a8] + HypertextLiteral v0.9.5
  [b5f81e59] + IOCapture v0.2.5
  [3587e190] + InverseFunctions v0.1.17
  [92d709cd] + IrrationalConstants v0.2.4
  [1019f520] + JLFzf v0.1.11
  [692b3bcd] + JLLWrappers v1.7.0
  [682c06a0] + JSON v0.21.4
  [b964fa9f] + LaTeXStrings v1.4.0
  [23fbe1c1] + Latexify v0.16.8
  [2ab3a3ac] + LogExpFunctions v0.3.28
  [e6f89c97] + LoggingExtras v1.1.0
  [6c6e2e6c] + MIMEs v1.1.0
  [1914dd2f] + MacroTools v0.5.16
  [739be429] + MbedTLS v1.1.9
  [442fdcdd] + Measures v0.3.2
  [e1d29d7a] + Missings v1.2.0
  [77ba4419] + NaNMath v1.0.3
  [4d8831e6] + OpenSSL v1.5.0
  [bac558e1] + OrderedCollections v1.8.1
  [90014a1f] + PDMats v0.11.31
  [69de0a69] + Parsers v2.8.3
  [ccf2f8ad] + PlotThemes v3.3.0
  [995b91a9] + PlotUtils v1.4.3
  [91a5bcdd] + Plots v1.40.14
  [7f904dfe] + PlutoUI v0.7.64
  [aea7be01] + PrecompileTools v1.2.1
  [21216c6a] + Preferences v1.4.3
  [43287f4e] + PtrArrays v1.3.0
  [1fd47b50] + QuadGK v2.11.2
  [3cdcf5f2] + RecipesBase v1.3.4
  [01d81517] + RecipesPipeline v0.6.12
  [189a3867] + Reexport v1.2.2
  [05181044] + RelocatableFolders v1.0.1
  [ae029012] + Requires v1.3.1
  [79098fc4] + Rmath v0.8.0
  [6c6a2e73] + Scratch v1.2.1
  [992d4aef] + Showoff v1.0.3
  [777ac1f9] + SimpleBufferStream v1.2.0
  [a2af1166] + SortingAlgorithms v1.2.1
  [276daf66] + SpecialFunctions v2.5.1
  [860ef19b] + StableRNGs v1.0.3
  [82ae8749] + StatsAPI v1.7.1
  [2913bbd2] + StatsBase v0.34.4
  [4c63d2b9] + StatsFuns v1.4.0
  [62fd8b95] + TensorCore v0.1.1
  [3bb67fe8] + TranscodingStreams v0.11.3
  [410a4b4d] + Tricks v0.1.10
  [5c2747f8] + URIs v1.5.2
  [1cfade01] + UnicodeFun v0.4.1
  [1986cc42] + Unitful v1.23.1
  [45397f5d] + UnitfulLatexify v1.7.0
  [41fe7b60] + Unzip v0.2.0
  [6e34b625] + Bzip2_jll v1.0.9+0
  [83423d85] + Cairo_jll v1.18.5+0
  [ee1fde0b] + Dbus_jll v1.16.2+0
  [2702e6a9] + EpollShim_jll v0.0.20230411+1
  [2e619515] + Expat_jll v2.6.5+0
  [b22a6f82] + FFMPEG_jll v4.4.4+1
  [a3f928ae] + Fontconfig_jll v2.16.0+0
  [d7e528f0] + FreeType2_jll v2.13.4+0
  [559328eb] + FriBidi_jll v1.0.17+0
  [0656b61e] + GLFW_jll v3.4.0+2
  [d2c73de3] + GR_jll v0.73.6+0
  [78b55507] + Gettext_jll v0.21.0+0
  [7746bdde] + Glib_jll v2.84.0+0
  [3b182d85] + Graphite2_jll v1.3.15+0
  [2e76f6c2] + HarfBuzz_jll v8.5.1+0
  [aacddb02] + JpegTurbo_jll v3.1.1+0
  [c1c5ebd0] + LAME_jll v3.100.2+0
  [88015f11] + LERC_jll v3.0.0+1
  [1d63c593] + LLVMOpenMP_jll v18.1.8+0
  [dd4b983a] + LZO_jll v2.10.3+0
  [e9f186c6] + Libffi_jll v3.4.7+0
  [7e76a0d4] + Libglvnd_jll v1.7.1+1
  [94ce4f54] + Libiconv_jll v1.18.0+0
  [4b2f31a3] + Libmount_jll v2.41.0+0
  [89763e89] + Libtiff_jll v4.5.1+1
  [38a345b3] + Libuuid_jll v2.41.0+0
  [e7412a2a] + Ogg_jll v1.3.5+1
  [458c3c95] + OpenSSL_jll v3.5.0+0
  [efe28fd5] + OpenSpecFun_jll v0.5.6+0
  [91d4177d] + Opus_jll v1.3.3+0
  [36c8627f] + Pango_jll v1.56.3+0
  [30392449] + Pixman_jll v0.44.2+0
  [c0090381] + Qt6Base_jll v6.7.1+1
  [f50d1b31] + Rmath_jll v0.5.1+0
  [a44049a8] + Vulkan_Loader_jll v1.3.243+0
  [a2964d1f] + Wayland_jll v1.23.1+0
  [2381bf8a] + Wayland_protocols_jll v1.44.0+0
  [02c8fc9c] + XML2_jll v2.13.6+1
  [ffd25f8a] + XZ_jll v5.8.1+0
  [f67eecfb] + Xorg_libICE_jll v1.1.2+0
  [c834827a] + Xorg_libSM_jll v1.2.6+0
  [4f6342f7] + Xorg_libX11_jll v1.8.12+0
  [0c0b7dd1] + Xorg_libXau_jll v1.0.13+0
  [935fb764] + Xorg_libXcursor_jll v1.2.4+0
  [a3789734] + Xorg_libXdmcp_jll v1.1.6+0
  [1082639a] + Xorg_libXext_jll v1.3.7+0
  [d091e8ba] + Xorg_libXfixes_jll v6.0.1+0
  [a51aa0fd] + Xorg_libXi_jll v1.8.3+0
  [d1454406] + Xorg_libXinerama_jll v1.1.6+0
  [ec84b674] + Xorg_libXrandr_jll v1.5.5+0
  [ea2f1a96] + Xorg_libXrender_jll v0.9.12+0
  [c7cfdc94] + Xorg_libxcb_jll v1.17.1+0
  [cc61e674] + Xorg_libxkbfile_jll v1.1.3+0
  [e920d4aa] + Xorg_xcb_util_cursor_jll v0.1.4+0
  [12413925] + Xorg_xcb_util_image_jll v0.4.1+0
  [2def613f] + Xorg_xcb_util_jll v0.4.1+0
  [975044d2] + Xorg_xcb_util_keysyms_jll v0.4.1+0
  [0d47668e] + Xorg_xcb_util_renderutil_jll v0.3.10+0
  [c22f9ab0] + Xorg_xcb_util_wm_jll v0.4.2+0
  [35661453] + Xorg_xkbcomp_jll v1.4.7+0
  [33bec58e] + Xorg_xkeyboard_config_jll v2.44.0+0
  [c5fb5394] + Xorg_xtrans_jll v1.6.0+0
  [3161d3a3] + Zstd_jll v1.5.7+1
  [35ca27e7] + eudev_jll v3.2.14+0
  [214eeab7] + fzf_jll v0.61.1+0
  [a4ae2306] + libaom_jll v3.11.0+0
  [0ac62f75] + libass_jll v0.15.2+0
  [1183f4f0] + libdecor_jll v0.2.2+0
  [2db6ffa8] + libevdev_jll v1.13.4+0
  [f638f0a6] + libfdk_aac_jll v2.0.3+0
  [36db933b] + libinput_jll v1.28.1+0
  [b53b4c65] + libpng_jll v1.6.49+0
  [f27f6e37] + libvorbis_jll v1.3.7+2
  [009596ad] + mtdev_jll v1.1.7+0
  [1270edf5] + x264_jll v2021.5.5+0
  [dfaa095f] + x265_jll v3.5.0+0
  [d8fb68d0] + xkbcommon_jll v1.8.1+0
  [0dad84c5] + ArgTools
  [56f22d72] + Artifacts
  [2a0f44e3] + Base64
  [ade2ca70] + Dates
  [8bb1440f] + DelimitedFiles
  [f43a241f] + Downloads
  [7b1f6079] + FileWatching
  [b77e0a4c] + InteractiveUtils
  [b27032c2] + LibCURL
  [76f85450] + LibGit2
  [8f399da3] + Libdl
  [37e2e46d] + LinearAlgebra
  [56ddb016] + Logging
  [d6f4376e] + Markdown
  [a63ad114] + Mmap
  [ca575930] + NetworkOptions
  [44cfe95a] + Pkg
  [de0858da] + Printf
  [3fa0cd96] + REPL
  [9a3f8284] + Random
  [ea8e919c] + SHA
  [9e88b42a] + Serialization
  [6462fe0b] + Sockets
  [2f01184e] + SparseArrays
  [10745b16] + Statistics
  [4607b0f0] + SuiteSparse
  [fa267f1f] + TOML
  [a4e569a6] + Tar
  [8dfed614] + Test
  [cf7118a7] + UUIDs
  [4ec0a83e] + Unicode
  [e66e0078] + CompilerSupportLibraries_jll
  [deac9b47] + LibCURL_jll
  [29816b5a] + LibSSH2_jll
  [c8ffd9c3] + MbedTLS_jll
  [14a3606d] + MozillaCACerts_jll
  [4536629a] + OpenBLAS_jll
  [05823500] + OpenLibm_jll
  [efcefdf7] + PCRE2_jll
  [83775a58] + Zlib_jll
  [8e850b90] + libblastrampoline_jll
  [8e850ede] + nghttp2_jll
  [3f19e933] + p7zip_jll
8.0 s
html"""
<iframe width="100%" height="300" src="https://www.youtube.com/embed/Gi4ZZVS2GLA" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>
"""
👀 Reading hidden code
120 μs

Before working on the homework, make sure that you have watched the first lecture on climate modeling 👆. We have included the important functions from this lecture notebook in the next cell. Feel free to have a look!

👀 Reading hidden code
259 μs
221.2
Model.A
👀 Reading hidden code
14.2 μs
Main.workspace#3.Model
module Model

const S = 1368; # solar insolation [W/m^2] (energy per unit time per unit area)
const α = 0.3; # albedo, or planetary reflectivity [unitless]
const B = -1.3; # climate feedback parameter [W/m^2/°C],
const T0 = 14.; # preindustrial temperature [°C]

absorbed_solar_radiation(; α=α, S=S) = S*(1 - α)/4; # [W/m^2]
outgoing_thermal_radiation(T; A=A, B=B) = A - B*T;

const A = S*(1. - α)/4 + B*T0; # [W/m^2].

greenhouse_effect(CO2; a=a, CO2_PI=CO2_PI) = a*log(CO2/CO2_PI);

const a = 5.0; # CO2 forcing coefficient [W/m^2]
const CO2_PI = 280.; # preindustrial CO2 concentration [parts per million; ppm];
CO2_const(t) = CO2_PI; # constant CO2 concentrations

const C = 51.; # atmosphere and upper-ocean heat capacity [J/m^2/°C]

function timestep!(ebm)
append!(ebm.T, ebm.T[end] + ebm.Δt*tendency(ebm));
append!(ebm.t, ebm.t[end] + ebm.Δt);
end;

tendency(ebm) = (1. /ebm.C) * (
+ absorbed_solar_radiation(α=ebm.α, S=ebm.S)
- outgoing_thermal_radiation(ebm.T[end], A=ebm.A, B=ebm.B)
+ greenhouse_effect(ebm.CO2(ebm.t[end]), a=ebm.a, CO2_PI=ebm.CO2_PI)
);

begin
mutable struct EBM
T::Array{Float64, 1}
t::Array{Float64, 1}
Δt::Float64
CO2::Function
C::Float64
a::Float64
A::Float64
B::Float64
CO2_PI::Float64
α::Float64
S::Float64
end;
# Make constant parameters optional kwargs
EBM(T::Array{Float64, 1}, t::Array{Float64, 1}, Δt::Real, CO2::Function;
C=C, a=a, A=A, B=B, CO2_PI=CO2_PI, α=α, S=S) = (
EBM(T, t, Δt, CO2, C, a, A, B, CO2_PI, α, S)
);
# Construct from float inputs for convenience
EBM(T0::Real, t0::Real, Δt::Real, CO2::Function;
C=C, a=a, A=A, B=B, CO2_PI=CO2_PI, α=α, S=S) = (
EBM(Float64[T0], Float64[t0], Δt, CO2;
C=C, a=a, A=A, B=B, CO2_PI=CO2_PI, α=α, S=S);
);
end;

begin
function run!(ebm::EBM, end_year::Real)
while ebm.t[end] < end_year
timestep!(ebm)
end
end;
run!(ebm) = run!(ebm, 200.) # run for 200 years by default
end




CO2_hist(t) = CO2_PI * (1 .+ fractional_increase(t));
fractional_increase(t) = ((t .- 1850.)/220).^3;

begin
CO2_RCP26(t) = CO2_PI * (1 .+ fractional_increase(t) .* min.(1., exp.(-((t .-1850.).-170)/100))) ;
RCP26 = EBM(T0, 1850., 1., CO2_RCP26)
run!(RCP26, 2100.)
CO2_RCP85(t) = CO2_PI * (1 .+ fractional_increase(t) .* max.(1., exp.(((t .-1850.).-170)/100)));
RCP85 = EBM(T0, 1850., 1., CO2_RCP85)
run!(RCP85, 2100.)
end

end
👀 Reading hidden code
179 ms

Exercise 1 - policy goals under uncertainty

A recent ground-breaking review paper produced the most comprehensive and up-to-date estimate of the climate feedback parameter, which they find to be

BN(1.3,0.4),

i.e. our knowledge of the real value is normally distributed with a mean value B=1.3 W/m²/K and a standard deviation σ=0.4 W/m²/K. These values are not very intuitive, so let us convert them into more policy-relevant numbers.

Definition: Equilibrium climate sensitivity (ECS) is defined as the amount of warming ΔT caused by a doubling of CO₂ (e.g. from the pre-industrial value 280 ppm to 560 ppm), at equilibrium.

At equilibrium, the energy balance model equation is:

0=S(1α)4(ABTeq)+aln(2COPICOPI)

From this, we subtract the preindustrial energy balance, which is given by:

0=S(1α)4(ABT0),

The result of this subtraction, after rearranging, is our definition of ECS:

ECSTeqT0=aln(2)B

👀 Reading hidden code
709 μs
ECS(; B=B̅, a=Model.a) = -a*log(2.)./B;
👀 Reading hidden code
1.6 ms

The plot below provides an example of an "abrupt 2xCO₂" experiment, a classic experimental treatment method in climate modelling which is used in practice to estimate ECS for a particular model. (Note: in complicated climate models the values of the parameters a and B are not specified a priori, but emerge as outputs of the simulation.)

The simulation begins at the preindustrial equilibrium, i.e. a temperature T0=14°C is in balance with the pre-industrial CO₂ concentration of 280 ppm until CO₂ is abruptly doubled from 280 ppm to 560 ppm. The climate responds by warming rapidly, and after a few hundred years approaches the equilibrium climate sensitivity value, by definition.

👀 Reading hidden code
360 μs
👀 Reading hidden code
3.9 s

B= -1.301

👀 Reading hidden code
372 ms
Exercise 1.1 - Develop understanding for feedbacks and climate sensitivity
👀 Reading hidden code
232 μs

👉 Change the value of B using the slider above. What does it mean for a climate system to have a more negative value of B? Explain why we call B the climate feedback parameter.

👀 Reading hidden code
278 μs

Hello world!

observations_from_changing_B = md"""
Hello world!
"""
👀 Reading hidden code
25.7 ms

👉 What happens when B is greater than or equal to zero?

👀 Reading hidden code
186 μs

Hello world!

observations_from_nonnegative_B = md"""
Hello world!
"""
👀 Reading hidden code
204 μs

Reveal answer:

👀 Reading hidden code
37.0 ms
👀 Reading hidden code
6.5 ms

👉 Create a graph to visualize ECS as a function of B.

👀 Reading hidden code
180 μs

👀 Reading hidden code
59.9 μs

Exercise 1.2 - Doubling CO₂

To compute ECS, we doubled the CO₂ in our atmosphere. This factor 2 is not entirely arbitrary: without substantial effort to reduce CO₂ emissions, we are expected to at least double the CO₂ in our atmosphere by 2100.

Right now, our CO₂ concentration is 415 ppm – 1.482 times the pre-industrial value of 280 ppm from 1850.

The CO₂ concentrations in the future depend on human action. There are several models for future concentrations, which are formed by assuming different policy scenarios. A baseline model is RCP8.5 - a "worst-case" high-emissions scenario. In our notebook, this model is given as a function of t.

👀 Reading hidden code
21.1 ms
plot(t, Model.CO2_RCP85.(t),
ylim=(0,1200), ylabel="CO2 concentration [ppm]")
👀 Reading hidden code
180 ms

👉 In what year are we expected to have doubled the CO₂ concentration, under policy scenario RCP8.5?

👀 Reading hidden code
187 μs
missing
expected_double_CO2_year = let
missing
end
👀 Reading hidden code
16.7 μs

Hint

The function findfirst might be helpful.

👀 Reading hidden code
45.5 ms

Exercise 1.3 - Uncertainty in B

The climate feedback parameter B is not something that we can control– it is an emergent property of the global climate system. Unfortunately, B is also difficult to quantify empirically (the relevant processes are difficult or impossible to observe directly), so there remains uncertainty as to its exact value.

A value of B close to zero means that an increase in CO₂ concentrations will have a larger impact on global warming, and that more action is needed to stay below a maximum temperature. In answering such policy-related question, we need to take the uncertainty in B into account. In this exercise, we will do so using a Monte Carlo simulation: we generate a sample of values for B, and use these values in our analysis.

👀 Reading hidden code
379 μs
0.4
B̅ = -1.3; σ = 0.4
👀 Reading hidden code
14.4 μs
👀 Reading hidden code
61.0 μs
histogram(B_samples, size=(600, 250), label=nothing, xlabel="B [W/m²/K]", ylabel="samples")
👀 Reading hidden code
2.1 s

👉 Generate a probability distribution for the ECS based on the probability distribution function for B above. Plot a histogram.

👀 Reading hidden code
191 μs
ECS_samples = ECS.(B=B_samples)
👀 Reading hidden code
45.3 ms
# ECS_samples = missing
👀 Reading hidden code
13.5 μs

Answer:

👀 Reading hidden code
217 μs
histogram(ECS_samples, xlims=(0, 8), size=(500, 240))
👀 Reading hidden code
79.1 ms

It looks like the ECS distribution is not normally distributed, even though B is.

👉 How does ECS(B) compare to ECS(B)? What is the probability that ECS(B) lies above ECS(B)?

👀 Reading hidden code
324 μs
2.6393298091640482
ecs_of_mean = ECS(B=mean(B_samples))
👀 Reading hidden code
24.7 μs
3.0136345332325356
mean_of_ecs = mean(ECS.(B=B_samples))
👀 Reading hidden code
2.7 ms
0.5011002200440088
sum(ECS_samples) do e
e > ecs_of_mean
end / length(ECS_samples)
👀 Reading hidden code
27.7 μs
0.34266853370674133
sum(ECS_samples) do e
e > mean_of_ecs
end / length(ECS_samples)
👀 Reading hidden code
25.8 μs

👉 Does accounting for uncertainty in feedbacks make our expectation of global warming better (less implied warming) or worse (more implied warming)?

👀 Reading hidden code
190 μs

Hello world!

observations_from_the_order_of_averaging = md"""
Hello world!
"""
👀 Reading hidden code
228 μs

Exercise 1.5 - Running the model

In the lecture notebook we introduced a mutable struct EBM (energy balance model), which contains:

  • the parameters of our climate simulation (C, a, A, B, CO2_PI, α, S, see details below)

  • a function CO2, which maps a time t to the concentrations at that year. For example, we use the function t -> 280 to simulate a model with concentrations fixed at 280 ppm.

EBM also contains the simulation results, in two arrays:

  • T is the array of tempartures (°C, Float64).

  • t is the array of timestamps (years, Float64), of the same size as T.

👀 Reading hidden code
781 μs

Properties of an EBM obect:

NameDescription
ALinearized outgoing thermal radiation: offset [W/m²]
BLinearized outgoing thermal radiation: slope. or: climate feedback parameter [W/m²/°C]
αPlanet albedo, 0.0-1.0 [unitless]
SSolar insulation [W/m²]
CAtmosphere and upper-ocean heat capacity [J/m²/°C]
aCO₂ forcing effect [W/m²]
CO2_PIPre-industrial CO₂ concentration [ppm]
👀 Reading hidden code
117 μs

You can set up an instance of EBM like so:

👀 Reading hidden code
189 μs
empty_ebm = Model.EBM(
14.0, # initial temperature
1850, # initial year
1, # Δt
t -> 280.0, # CO2 function
)
👀 Reading hidden code
11.1 ms

Have look inside this object. We see that T and t are initialized to a 1-element array.

Let's run our model:

👀 Reading hidden code
233 μs
simulated_model = let
ebm = Model.EBM(14.0, 1850, 1, t -> 280.0)
Model.run!(ebm, 2020)
ebm
end
👀 Reading hidden code
121 μs

Again, look inside simulated_model and notice that T and t have accumulated the simulation results.

In this simulation, we used T0 = 14 and CO2 = t -> 280, which is why T is constant during our simulation. These parameters are the default, pre-industrial values, and our model is based on this equilibrium.

👉 Run a simulation with policy scenario RCP8.5, and plot the computed temperature graph. What is the global temperature at 2100?

👀 Reading hidden code
314 μs
missing
simulated_rcp85_model = let
missing
end
👀 Reading hidden code
18.2 μs

👀 Reading hidden code
59.3 μs

Additional parameters can be set using keyword arguments. For example:

Model.EBM(14, 1850, 1, t -> 280.0; B=-2.0)

Creates the same model as before, but with B = -2.0.

👀 Reading hidden code
260 μs

👉 Write a function temperature_response that takes a function CO2 and an optional value B as parameters, and returns the temperature at 2100 according to our model.

👀 Reading hidden code
208 μs
temperature_response (generic function with 2 methods)
function temperature_response(CO2::Function, B::Float64=-1.3)
return missing
end
👀 Reading hidden code
605 μs
missing
temperature_response(t -> 280)
👀 Reading hidden code
15.2 μs
missing
temperature_response(Model.CO2_RCP85)
👀 Reading hidden code
17.8 μs
missing
temperature_response(Model.CO2_RCP85, -1.0)
👀 Reading hidden code
17.5 μs

Exercise 1.6 - Application to policy relevant questions

We talked about two emissions scenarios: RCP2.6 (strong mitigation - controlled CO2 concentrations) and RCP8.5 (no mitigation - high CO2 concentrations). These are given by the following functions:

👀 Reading hidden code
327 μs
Model.CO2_RCP26(t_scenario_test), Model.CO2_RCP85(t_scenario_test)
👀 Reading hidden code
20.3 ms
1850
@bind t_scenario_test Slider(t; show_value=true, default=1850)
👀 Reading hidden code
104 ms
1850:2100
👀 Reading hidden code
12.3 μs

We are interested in how the uncertainty in our input B (the climate feedback paramter) propagates through our model to determine the uncertainty in our output T(t), for a given emissions scenario. The goal of this exercise is to answer the following by using Monte Carlo Simulation for uncertainty propagation:

👉 What is the probability that we see more than 2°C of warming by 2100 under the low-emissions scenario RCP2.6? What about under the high-emissions scenario RCP8.5?

👀 Reading hidden code
419 μs

👀 Reading hidden code
59.6 μs

Exercise 2 - How did Snowball Earth melt?

In lecture 21 (see below), we discovered that increases in the brightness of the Sun are not sufficient to explain how Snowball Earth eventually melted.

👀 Reading hidden code
322 μs
👀 Reading hidden code
112 μs

We talked about a second theory – a large increase in CO₂ (by volcanoes) could have caused a strong enough greenhouse effect to melt the Snowball. If we imagine that the CO₂ then decreased (e.g. by getting sequestered by the now liquid ocean), we might be able to explain how we transitioned from a hostile Snowball Earth to today's habitable "Waterball" Earth.

In this exercise, you will estimate how much CO₂ would be needed to melt the Snowball and visualize a possible trajectory for Earth's climate over the past 700 million years by making an interactive bifurcation diagram.

Exercise 2.1

In the lecture notebook (video above), we had a bifurcation diagram of S (solar insolation) vs T (temperature). We increased S, watched our point move right in the diagram until we found the tipping point. This time we will do the same, but we vary the CO₂ concentration, and keep S fixed at its default (present day) value.

👀 Reading hidden code
479 μs

Below we have an empty diagram, which is already set up with a CO₂ vs T diagram, with a logirthmic horizontal axis. Now it's your turn! We have written some pointers below to help you, but feel free to do it your own way.

👀 Reading hidden code
223 μs
let
p = plot(
xlims=(CO2min, CO2max), ylims=(-55, 75),
xaxis=:log,
xlabel="CO2 concentration [ppm]",
ylabel="Global temperature T [°C]",
title="Earth's CO2 concentration bifurcation diagram",
legend=:topleft
)
add_cold_hot_areas!(p)
add_reference_points!(p)
# your code here
# TODO:
step_model!(ebm, CO2)
plot!(p,
[ebm.CO2(ebm.t[end])], [ebm.T[end]],
label=nothing,
color=:black,
shape=:circle,
)
end |> as_svg
👀 Reading hidden code
991 ms

We used two helper functions:

👀 Reading hidden code
179 μs
add_cold_hot_areas! (generic function with 1 method)
👀 Reading hidden code
1.3 ms
add_reference_points! (generic function with 1 method)
👀 Reading hidden code
1.1 ms

👉 Create a slider for CO2 between CO2min and CO2max. Just like the horizontal axis of our plot, we want the slider to be logarithmic.

👀 Reading hidden code
243 μs
281.8382931264455
CO2 = 10^log_CO2
👀 Reading hidden code
15.1 μs
begin
@bind log_CO2 Slider(log10(CO2min):0.01:log10(CO2max); default=log10(Model.CO2_PI))
end
👀 Reading hidden code
19.2 ms

Hint

@bind log_CO2 Slider(❓)
CO2 = 10^log_CO2
👀 Reading hidden code
184 μs
ebm = Model.EBM(Tneo, 0., 5., Model.CO2_const)
👀 Reading hidden code
9.4 ms

👉 Write a function step_model! that takes an existing ebm and new_CO2, which performs a step of our interactive process:

  • Reset the model by setting the ebm.t and ebm.T arrays to a single element. Which value?

  • Assign a new function to ebm.CO2. What function?

  • Run the model.

👀 Reading hidden code
531 μs
step_model! (generic function with 1 method)
function step_model!(ebm::Model.EBM, new_CO2::Real)
ebm.T = [ebm.T[end]]
ebm.t = [0]
ebm.CO2 = t -> new_CO2
Model.run!(ebm, 500)
ebm
end
👀 Reading hidden code
1.2 ms
# function step_model!(ebm::Model.EBM, CO2::Real)
# # your code here
# return ebm
# end
👀 Reading hidden code
13.8 μs

👉 Inside the plot cell, call the function step_model!.

👀 Reading hidden code
196 μs
Parameters
👀 Reading hidden code
176 μs
10
👀 Reading hidden code
11.4 μs
1000000
👀 Reading hidden code
10.8 μs
-48
👀 Reading hidden code
10.0 μs

The albedo feedback is implemented by the methods below:

👀 Reading hidden code
200 μs
α (generic function with 1 method)
function α(T; α0=Model.α, αi=0.5, ΔT=10.)
if T < -ΔT
return αi
elseif -ΔT <= T < ΔT
return αi + (α0-αi)*(T+ΔT)/(2ΔT)
elseif T >= ΔT
return α0
end
end
👀 Reading hidden code
2.0 ms
function Model.timestep!(ebm)
ebm.α = α(ebm.T[end]) # Added this line
append!(ebm.T, ebm.T[end] + ebm.Δt*Model.tendency(ebm));
append!(ebm.t, ebm.t[end] + ebm.Δt);
end
👀 Reading hidden code
875 μs

If you like, make the visualization more informative! Like in the lecture notebook, you could add a trail behind the black dot, or you could plot the stable and unstable branches. It's up to you!

👀 Reading hidden code
189 μs

Exercise 2.2

👉 Find the lowest CO₂ concentration necessary to melt the Snowball, programatically.

👀 Reading hidden code
260 μs
missing
co2_to_melt_snowball = let
missing
end
👀 Reading hidden code
16.5 μs

Hint

Use a condition on the albedo or temperature to check whether the Snowball has melted.

👀 Reading hidden code
260 μs

Hint

Start by writing a function equilibrium_temperature(CO2) which creates a new EBM at the Snowball Earth temperature T = -48 and returns the final temperature for a given CO2 level.

👀 Reading hidden code
1.6 ms

Exercise XX: Lecture transcript

(MIT students only)

Please see the link for hw 9 transcript document on Canvas. We want each of you to correct about 500 lines, but don’t spend more than 20 minutes on it. See the the beginning of the document for more instructions. :point_right: Please mention the name of the video(s) and the line ranges you edited:

👀 Reading hidden code
422 μs

Abstraction, lines 1-219; Array Basics, lines 1-137; Course Intro, lines 1-144 (for example)

lines_i_edited = md"""
Abstraction, lines 1-219; Array Basics, lines 1-137; Course Intro, lines 1-144 (_for example_)
"""
👀 Reading hidden code
297 μs

Before you submit

Remember to fill in your name and Kerberos ID at the top of this notebook.

👀 Reading hidden code
852 μs

Function library

Just some helper functions used in the notebook.

👀 Reading hidden code
246 μs
hint (generic function with 1 method)
👀 Reading hidden code
491 μs
almost (generic function with 1 method)
👀 Reading hidden code
457 μs
still_missing (generic function with 2 methods)
👀 Reading hidden code
813 μs
keep_working (generic function with 2 methods)
👀 Reading hidden code
809 μs
👀 Reading hidden code
10.1 ms
correct (generic function with 2 methods)
👀 Reading hidden code
725 μs
not_defined (generic function with 1 method)
👀 Reading hidden code
805 μs
TODO
👀 Reading hidden code
33.4 μs